This IPython notebook builds on the previous blog post which described how to simulate and plot the result of a shaped pulse on magnetization in an NMR experiment. For the purpose of teaching1, it can also be useful to visualize the effect of this pulse on the magnetization at each discrete time step of the pulse.
Visualizing the time-dependent course of magnetization requires the javascript viewer developed by Jake Vanderplas described here and further demonstrated here. This library has to be installed somewhere in the path searched by Python. NymPy and Matplotlib are also required.
from JSAnimation.IPython_display import display_animation
from matplotlib.animation import FuncAnimation
The steps required to created the Reburp shaped pulse and calculate its propagator have been described in a previous post and are available in the downloadable and static version of this notebook linked to below. Thus, I have jumped directly to the data visualization.
xyzdata = pulseprop_pulse_steps(relativeomega, pulseShapeInten, pulseShapePhase, gammaB1max, nu1maxdt, vectorComponent, n_pulse, n_freq)
The entire range of frequencies and pulse steps isn't needed for the animation itself (the entire frequency range is needed for the first graph below, however), so they can be reduced to a smaller set at regular intervals.
freq_skip = 50
pulse_skip = 10
xyzdata_S = xyzdata[:,::freq_skip,::pulse_skip]
relativeomega_S = relativeomega[::freq_skip]
A colormap can help visualize a range of data plotted simultaneously in Matplotlib. For the animation, it would be useful to use this colormap as a way to indicate which frequencies are at the limit or outside of the bandwidth of the shaped pulse. This can be accomplished using a custom colormap which is then mapped to the resonance offset range.
from matplotlib.colors import LinearSegmentedColormap, Normalize
from matplotlib.cm import ScalarMappable
def omega_colors(relativeomega):
# First create the colormap where extremes are red, middle is blue
cdict = {'red': ((0.0, 1.0, 1.0),
(0.35, 0.0, 0.0),
(0.65, 0.0, 0.0),
(1.0, 1.0, 1.0)),
'green': ((0.0, 0.0, 0.0),
(1.0, 0.0, 0.0)),
'blue': ((0.0, 0.0, 0.0),
(0.35, 1.0, 1.0),
(0.65, 1.0, 1.0),
(1.0, 0.0, 0.0)) }
rbr = LinearSegmentedColormap('RedBlueRed', cdict)
# Normalize the colormap to the extremes of the frequency range
cnorm = Normalize(vmin=np.min(relativeomega),vmax=np.max(relativeomega))
# Generate a blended map over the entire frequency range
# Then get colors for each offset value
scalarMap = ScalarMappable(norm=cnorm, cmap=rbr)
colorList = [scalarMap.to_rgba(x) for x in relativeomega]
return colorList, scalarMap
Plot the final magnetization and also visualize the frequencies that will be animated.
fig = plt.figure()
fig.set_facecolor((1.0,1.0,1.0,1.0))
ax = plt.subplot(111)
ax.plot(relativeomega,xyzdata[0,:,-1],color=(0.8,0.8,0.8,1.0),\
label='Mx',linewidth=3.0)
ax.plot(relativeomega,xyzdata[1,:,-1],color=(0.4,0.4,0.4,1.0),\
label='My',linewidth=3.0)
ax.plot(relativeomega,xyzdata[2,:,-1],color=(0.0,0.0,0.0,1.0),\
label='Mz',linewidth=3.0)
ax.set_xlabel('Omega (Hz)')
ax.set_ylabel('Relative Intensity')
ax.set_xlim(1.1*offset)
ax.set_ylim(-1.1,1.1)
ax.legend()
# Generate the color list using the above function to visualize the frequencies used
# The scalarMap will be used in the animation plot
colorList, scalarMap = omega_colors(relativeomega_S)
for x in range(len(relativeomega_S)):
ax.plot([relativeomega_S[x],relativeomega_S[x]],[-1.1,1.1],color=colorList[x])
Matplotlib doesn't have a specific type of 3D axis for plotting vectors on a sphere, so I created a function that will setup the 3D axis.
from mpl_toolkits.mplot3d import Axes3D, proj3d
from mpl_toolkits.mplot3d import proj3d
def circular_axes(fig, title):
# Create a 3D figure
ax = fig.gca(projection='3d')
# Draw axes lines
ax.plot([-1,1],[0,0],[0,0],'gray')
ax.plot([0,0],[-1,1],[0,0],'gray')
ax.plot([0,0],[0,0],[-1,1],'gray')
# Draw a gray circle in xy-plane
npoints = 100
ang = np.linspace(0,2*np.pi,npoints)
xcir = np.cos(ang)
ycir = np.sin(ang)
zcir = np.zeros((npoints,))
ax.plot(xcir,ycir,zcir,'gray')
# Uncomment to draw circular axes limits for other planes
#ax.plot(zcir,xcir,ycir,'gray') # Circular permutations!
#ax.plot(ycir,zcir,xcir,'gray')
# Hide the gray background and set perspective
ax.set_axis_off()
ax.view_init(20,30)
# Can't use axis_label commands because they will label in the wrong spot
# Annotation commands only work on 2D plots, but can use in 3D plots if
# the 2D projection is calculated
xs_x, ys_x, _ = proj3d.proj_transform(1,0,0, ax.get_proj())
xs_y, ys_y, _ = proj3d.proj_transform(0,1,0, ax.get_proj())
xs_z, ys_z, _ = proj3d.proj_transform(0,0,1, ax.get_proj())
xs_t, ys_t, _ = proj3d.proj_transform(0.5,-1.0,1, ax.get_proj())
# Label the axes with the calculated 2D points
ax.annotate('X',xy=(xs_x,ys_x),xycoords='data',xytext=(-12,-10),\
textcoords='offset points',fontweight='bold')
ax.annotate('Y',xy=(xs_y,ys_y),xycoords='data',xytext=(10,-10),\
textcoords='offset points',fontweight='bold')
ax.annotate('Z',xy=(xs_z,ys_z),xycoords='data',xytext=(8,8),\
textcoords='offset points',fontweight='bold')
ax.annotate(title,xy=(xs_t,ys_t),xycoords='data',xytext=(-10,8),\
textcoords='offset points',fontweight='bold')
return ax
For vectors plotted on a sphere, I prefer the aesthetic of arrows to simple lines. This is an adaptation of a class written to plot arrows on three dimensional axes. I added a subroutine for updating the position of an existing arrow as well as a way to create a legend from the arrows if desired.
from matplotlib.patches import FancyArrowPatch
class Arrow3D(FancyArrowPatch):
def __init__(self, xs, ys, zs, *args, **kwargs):
FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
self._verts3d = xs, ys, zs
# Create a dummy line at the axes origin so the legend will work
ax = plt.gca() # Risky, but get_axes() doesn't work b/c axis not set yet
ax.plot3D([0,0],[0,0],[0,0],linewidth=1.0,\
label=kwargs['label'],color=kwargs['color'])
def draw(self, renderer):
xs3d, ys3d, zs3d = self._verts3d
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
FancyArrowPatch.draw(self, renderer)
# This will update the arrow position
def set_data3d(self, xs3d, ys3d, zs3d):
ax = self.get_axes()
rendererM = ax.get_proj()
xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, rendererM)
self._verts3d = xs3d, ys3d, zs3d
self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
And that's all the necessary setup required for the animation.
from matplotlib.colorbar import Colorbar
fig = plt.figure()
fig.set_facecolor((1.0,1.0,1.0,1.0))
# Create the fancy spherical axis
ax = circular_axes(fig,'Reburp Pulse')
# Set the starting position of all the arrows
arrowList = []
for x in range(len(relativeomega_S)):
arrow = Arrow3D([0,0],[0,0],[0,1],mutation_scale=20, lw=1.0, \
arrowstyle='-|>',color=colorList[x],label='%.0f Hz'%relativeomega_S[x])
arrowList.append(arrow)
ax.add_artist(arrow)
# Option 1: show a colorbar
scalarMap._A = []
cbar = fig.colorbar(scalarMap, shrink=.6, aspect=10)
cbar.ax.set_yticklabels(['%5.0f Hz'%x for x in relativeomega_S])
# Option 2: show a simple legend
# ax.legend()
# Function to update the arrow position for each pulse step during animation
def update_pos(time, data, arrowList):
xyz = data[:,:,time]
for x in range(len(arrowList)):
arrowList[x].set_data3d([0,xyz[0,x]],[0,xyz[1,x]],[0,xyz[2,x]])
# Pad the data with starting and finishing values for aesthetics
pad = 20
xyzdata_beg = np.tile(xyzdata_S[:,:,0][:,:,np.newaxis],(1,1,pad))
xyzdata_end = np.tile(xyzdata_S[:,:,-1][:,:,np.newaxis],(1,1,pad))
xyzdata_pad = np.append(xyzdata_beg,xyzdata_S,axis=-1)
xyzdata_pad = np.append(xyzdata_pad,xyzdata_end,axis=-1)
# Do the animation!
anim = FuncAnimation(fig,update_pos,frames=xyzdata_pad.shape[2], \
fargs=(xyzdata_pad, arrowList),interval=25)
display_animation(anim,default_mode='once')
Note that the arrows representing frequencies within the appropriate bandwidth are all inverted by the pulse. However, the arrows that reside on the edge of the bandwidth or outside of it are not inverted at all. This demonstrates the importance of careful calibration of shaped pulses and the importance of using gradients to clean up stray magnetization whenever possible.
I also think this is a cool animation that would be fun to use for teaching!
This post was written in an IPython notebook, which can be downloaded here, or viewed statically here.
-
And providing an excuse for me to play with animations in IPython notebooks. ↩